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PREFACE 


This  study  was  conducted  under  ILIR  Project  4A161 101A91D,  “Development  of  A 
New  Technique  for  Texture  Analysis.” 

The  work  described  in  this  research  note  represents  an  application  of  a  third-order 
co-occurrence  approach  toward  texture  analysis.  The  task  was  performed  under  the  super¬ 
vision  of  Dr.  Frederick  W.  Rohde,  Team  Leader,  Center  for  Physical  Sciences;  and  Dr. 
Robert  D.  Leighty,  Director,  Research  Institute. 

COL  Alan  L.  Laubscher,  CE,  was  the  Commander  and  Director  and  Mr.  Walter  E. 
Boge  was  the  Technical  Director  of  the  Engineer  Topographic  Laboratories  during  the 
report  preparation. 
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THIRD-ORDER  CO-OCCURRENCE  TEXTURE  ANALYSIS  APPLIED 
TO  SAMPLES  OF  HIGH  RESOLUTION  RADAR  IMAGERY 


INTRODUCTION 

The  purpose  of  this  research  note  is  to  show  the  application  of  a  third-order 
co-occurrence  texture  analysis  to  samples  of  high  resolution  synthetic  aperture  radar 
imagery.  In  the  past,  various  types  of  texture  analysis  techniques  have  been  applied  to 
digital  images.  One  of  the  most  popular  techniques  has  been  the  computation  of  the  gray 
level  co-occurrence  matrix.  Each  element  of  this  matrix  represents  the  number  of  times 
that  gray  level  £  occurs  next  to  gray  level  m  for  a  given  pixel  spacing  and  direction.  The 
concept  employed  in  this  report  is  a  natural  extension  of  this  technique.  In  a  third-order 
co-occurrence  analysis,  calculations  are  made  for  the  number  of  times  gray  level  £  occurs 
next  to  gray  level  m,  which  in  turn  occurs  next  to  gray  level  n.  An  interesting  aspect  of  this 
jump  from  second-  to  third-order  co-occurrence  is  that  now  all  three  pixels  do  not 
necessarily  have  to  lie  in  a  straight  line.  The  pixel  spacing  and  direction  between  gray 
level  £  and  gray  level  m  can  be  quite  different  from  the  spacing  and  direction  between 
gray  level  m  and  gray  level  n.  The  following  sections  will  present  a  discussion  of  the 
third-order  co-occurrence  technique  along  with  the  feature  selection  and  classification 
methods  applied  to  certain  samples  of  radar  imagery. 


METHODOLOGY 

This  section  develops  the  third-order  co-occurrence  method  and  discusses  briefly 
the  methods  used  for  feature  selection  and  classification.  The  third-order  approach  com¬ 
putes  the  number  of  times  in  a  given  image  that  gray  level  £  occurs  next  to  gray  level  m, 
which  in  turn  occurs  next  to  gray  level  n.  The  pixel  spacing  between  gray  level  £  and  gray 
1;>  m  ■>  tv  •.  >\v  Ax |  and  Ayj,  and  the  spacing  between  gray  level  £  and  gray 

level  n  will  be  denoted  by  Ax2  and  Ay2.  The  functional  notation  for  the  third-order 
co-occurrence  is 

N(  i,  j,  k|Ax,  ,  Ay,  ,  Ax2  ,  Ay2  ) 
where  i=£+l  j=m+l  k=n+l 


The  above  expression  will  be  shortened  to  Njj^,  where  the  values  for  the  four  spacing 
parameters  are  assumed.  Njjj,  can  be  visualized  as  a  third-order  array,  the  dimensions  of 
which  will  be  determined  by  the  total  number  of  gray  levels  in  the  digital  image.  If  the 
number  of  gray  levels  in  the  image  is  16,  then  Njj^  is  an  array  whose  dimensions  are  16  by 
16  by  16.  Various  measures  or  features  can  then  be  computed  from  Njjk  to  form  a  feature 
vector  that  can  be  used  for  classification  purposes.  The  13  features  given  below  are  cal¬ 
culated  from  Njjk  and  used  to  form  the  components  of  a  feature  vector  X. 
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A  general  expression  for  the  summation  of  terms  along  specific  diagonals  is  given  by  A 
below: 

Aab  “  ZZZ  Nijk 

i  j  k  J 
|i-j|-a  I i-k I -b 


Using  the  above  definition,  the  following  three  features  can  be  computed. 


A  -  Z  Z  Z  Nijk  -  x 

00  1  j  k  J  11 

|i-J 1-0  li-kl-0 


A01  -  2  Z  Z  HlJk  .  x 

01  1  j  k  12 

| i-k | -1 
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In  all  the  above  definitions,  the  summations  in  i,  j,  and  k  are  over  the  entire  range  of 
gray  levels  in  the  image.  If  there  are  16  gray  levels  in  the  image,  then  i,  j,  and  k  all  vary  from 
1  to  16.  There  is  nothing  particularly  unique  about  the  13  features  chosen  above.  Other 
features  could  also  be  computed  from  Njjk  and  perhaps  would  even  work  better  than  the 
set  chosen.  Once  the  form  of  the  feature  vector  X  has  been  set,  a  method  for  performing 
feature  selection  needs  to  be  considered. 

In  this  study  the  feature  selection  technique  came  from  the  fitid  of  discriminant 
analysis  of  statistics.  The  feature  selection  method  was  used  to  reduce  the  dimensionality  of 
X  from  13  to  2  and  at  the  same  time  optimize  the  separation  among  data  belonging  to 
different  classes.  The  feature  selection  operation  can  be  performed  by  using  a  linear 
transformation  as  follows: 


Y  =  AX 


where  X  is  the  original  feature  vector  with  dimensionality  13x1;  A  is  the  transformation 
matrix  of  dimensionality  2x13;  and  Y  is  the  transformed  feature  vector  with  the  dimen¬ 
sionality  2x1.  In  order  to  calculate  the  matrix  A,  use  was  made  of  the  within-class  and 
between-class  scatter  matrices.  A  within-class  scatter  matrix  shows  the  scatter  of  samples 
around  their  class  expected  vector  and  can  be  expressed  by 
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where  S,„  is  the  within-class  scatter  matrix,  P(w:)  is  the  a  priori  probability  of  the  i1*1  class, 
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w  th  ' 

C;  is  the  covariance  matrix  of  the  iUI  class,  and  Nt  is  the  total  number  of  classes.  A  between- 
class  scatter  matrix  can  be  defined  in  many  ways,  however,  the  following  definition  was  the 
one  utilized  here: 


( 1 6) 


where  is  the  between-class  scatter  matrix;  Cj  is  the  covariance  matrix  for  class  1;  C2 
is  the  covariance  matrix  for  class  2;  M  j  is  the  mean  vector  for  class  1 ;  M2  is  the  mean  vector 
for  class  2;  and  T  means  transpose.  This  definition  of  the  between-class  scatter  matrix  is 
valid  only  for  the  case  when  Nt  (the  total  number  of  classes)  is  equal  to  two.  In  order  to 
have  criteria  for  class  separability,  a  number  must  be  derived  from  the  scatter  matrices.  This 
number  should  increase  when  the  distances  between  points  belonging  to  different  classes  are 
increasing  or  when  the  distances  between  points  belonging  to  the  same  class  are  decreasing. 
One  criterion  is  the  use  of  Jj,  which  is  defined  as  follows: 
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The  feature  selection  problem  now  requires  that  we  select  the  particular  transformation 
matrix  A,  which  maximizes  the  value  of  J..  Fukunaga1  shows  that  A  is  made  up  of  the 
normalized  eigenvectors  of  the  matrix  S'1  Su . 

W  D 


A 


T 


(18) 


where  0 j  is  the  eigenvector  associated  with  the  largest  eigenvalue,  and  <£2  is  the  eigen¬ 
vector  associated  with  the  second  largest  eigenvalue.  Once  the  matrix  A  is  computed  from 
(18),  the  new  feature  vector  Y  can  be  computed  for  each  point  in  each  class.  After  feature 
selection  has  been  performed,  it  is  necessary  to  choose  a  classification  procedure. 

The  classification  procedure  used  in  this  work  was  the  pseudoinverse  technique.  This 
technique  is  a  linear  classifier,  which  attempts  to  solve  for  a  vector  ji  such  that 

if  _a^y  >  0  then  decide  class  1 

or 

if  a^y  <  0  then  decide  class  2 


where  the  vector  y  is  an  augmented  form  of  Y  as  shown  below : 


where  y,  and  y2  are  the  components  of  Y.  The  above  decision  function  is  valid  for  the  two 
class  problem.  It  can  easily  be  seen  that  the  vector  .a  for  our  case  has  three  components. 


'  Keinosukc  Fukunaga,  Introduction  to  Statistical  Pattern  Recognition ,  Academic  Press. 
1972. 
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A  solution  for  _a  can  be  found  by  forming  a  matrix  H  from  all  the  training  samples  taken 
from  the  two  classes.  Each  row  of  H  will  consist  of  a  sample  of  yT.  The  matrix  H  will 
therefore  have  the  following  form: 


1 

yn 

y2i 

1 

yn 

y22 

1 

yn 

y23 

1 

y  1  n 

y2n 

- 1 

"  y  1  n+ 1 

-y2n+l 

- 1 

' y 1 n+2 

■y2n+2 

- 1 

• 

• 

*  y 1 n+3 

• 

• 

’ y2n+3 

• 

• 

- 1 

' y  In+m 

• 

" y2n+m 

The  first  n  rows  of  H  consist  of  data  that  come  from  class  1 .  The  next  m  rows  of  H 
come  from  class  2  and  have  all  been  multiplied  by  -1.  The  second  subscript  on  the  y  values 
refers  to  the  sample  number.  An  equation  involving  H  and  £  is  given  below : 

Ha  =  b 

The  vector  has  n+m  components  all  of  which  are  arbitrarily  specified  positive  constants. 
Duda  and  Hart2  show  how  the  following  solution  for£  comes  from  equation  (20): 

a  =  (hTh)*1  HTb 

Once  the  three  components  of  a_  have  been  calculated  for  each  possible  pair  of  classes, 
then  the  classifier  has  been  completed  and  can  be  tested  against  the  original  training  set 
data  and  some  unknown  samples. 


2 

Richard  O.  Duda  and  Peter  E.  Hart,  Pattern  Classification  and  Scene  Analysis,  John  Wiley 
and  Sons,  1973. 
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INVESTIGATION 


The  third-order  co-occurrence  texture  calculations  were  applied  to  samples  of  high 
resolution  synthetic  aperture  radar  imagery  taken  over  the  Elizabeth  City,  North  Carolina 
area  with  the  UPD-4  radar  system.  Sections  of  the  radar  imagery  were  digitized  and  stored 
on  a  digital  disk  unit.  A  Lexidata  system  3400  display  processor  was  used  to  display  the 
images  on  a  cathode  ray  tube  and  to  take  700  samples  from  four  terrain  classes.  Each 
sample  consisted  of  a  32  by  32  pixel  element  window  located  within  a  homogeneous  section 
of  one  particular  terrain  class.  The  four  classes  considered  were  forests,  fields,  cities,  and 
water.  Of  the  700  samples  taken,  200  came  from  forests,  200  from  fields,  200  from  water, 
and  100  from  cities.  A  training  set  consisting  of  100  samples  from  each  of  the  four  classes 
was  used  to  compute  the  A  matrices  and  the  a.  vectors  for  each  of  the  six  possible  pairs  of 
classes.  Before  the  matrix  A  and  the  vector _a  could  be  calculated,  however,  it  was  necessary 
to  compute  Jj  as  a  function  of  several  values  of  Axj ,  Ayj ,  Ax2,  and  Ay2  in  an  attempt  to 
find  the  spacings  that  yield  a  significantly  large  value  for  Jj .  A  computer  program  was 
written  for  the  Hewlett-Packard  1000  computer  to  calculate  J,  as  a  function  of  the  spacing 
values.  This  program  was  also  used  to  compute  the  transformation  matrix  A.  A  listing  of 
this  program  along  with  the  subroutine  used  to  compute  the  third-order  co-occurrence 
array  is  given  in  appendix  A.  A  second  computer  program  was  written  to  calculate  the  new 
feature  vector  Y  and  to  compute  the  a. vector  from  the  components  of  Y.  A  listing  of  this 
second  program  is  provided  in  appendix  B.  A  third  computer  program  was  written  that  takes 
image  samples  and  makes  a  classification,  using  as  input  the  A  matrix  and  the  a^ vector  data 
calculated  in  the  previous  two  programs.  The  listing  for  this  third  program  is  given  in 
appendix  C. 


RESULTS 

In  this  section  some  results  of  numerical  calculations  are  presented.  Table  1  shows  the 
results  of  computing  the  value  of  Jj  for  each  of  the  six  possible  pairs  of  classes.  The  values 
for  the  spacing  parameters  (Axj,  Ayj,  Ax2,  Ay2)  shown  in  table  1  are  those  that  yielded  a 
maximum  value  for  J  j  for  each  pair  of  classes. 


Table  1.  Values  of  J  j  for  Third-Order  Features 


A*] 

Ayj 

Ax2 

Ay2 

Jl 

Forests  and  Fields 

1 

0 

2 

0 

49.78 

Forests  and  Cities 

2 

0 

2 

1 

23.05 

Cities  and  Fields 

1 

0 

2 

0 

169.21 

Forests  and  Water 

1 

0 

2 

0 

561.63 

Fields  and  Water 

1 

0 

2 

0 

196.66 

Cities  and  Water 

1 

0 

2 

0 

2790.16 

i. 


It  can  easily  be  seen  that  the  largest  three  values  of  J  t  all  occur  with  water  and  some 
other  class.  These  three  values  of  Jj  are  large  enough  to  say  that  the  data  from  these 
respective  classes  are  well  separated  and  clustered.  The  smallest  value  of  J  j  occurs  for  forests 
and  cities  and  indicates  that  the  classes  are  not  completely  separated.  Using  the  spacing 
values  given  in  table  1,  the  matrix  A  and  vector  a  were  computed  for  each  of  the  six  pairs 
of  classes,  as  shown  below: 


1.  Forests  and  Fields 


A  = 


-2.3x1  O'3  -0.996  -l.OlxlO'4  2.64x1  O'3  9.45x1  O'5  3.86x1  O'4  1.2x1 0‘5 
3.94x1 0‘3  0.999  2.78xl0‘S  -l.OlxlO'2  6.74xl0*3  I.68xl0'3  1.33xl0'6 


4. 12x10* 2  -3.89x1  O'5  -7.95xl0'2  -2.07xl0*4  -5.08xl04  -4.74xl0'4 
-4.43xl0‘2  -8.15xl0'6  -l.OOxlO'2  -9.40xl0'S  -7.7xl0'5  -2.28xl0'5 


V 


h 


-5.63792 

-19.91341 

-0.00121 


2.  Forests  and  Cities 

1.57xl0'2  -0.993  3.5xl0'5  -1.43x1  O'3  -2.07xI0"4  2.84xl0'3  2.21xl0'6 
2.4x1  O'3  -0.992  -4.34x1  O'5  1.07x1  O'2  1.35x1  O'3  -9.54x1  O'3  -lxl O'5 

-0.117  2.57xl0"4  -1.07xl0*2  -5.79xl0'5  1.39xl0'4  -8.72xl0"5 
0.101  —  1 .57x1 0*4  7.33x10' 2  4.36xI0'4  2.64xl0'4  4.31xl0'4 


a  = 


1.40889  “ 
113.96356 
-0.00306 
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3.  Cities  and  Fields 


-4.93x1  O’2  -1.27X10'1  3.6xlO'4  5.68xlO'3  -6.98xl0'4  7.44xl0'3  2.53xlO'5 
-5.1 1x1  O'3  9.98X10'1  2.1 1x10' 3  -1.46xl0"2  1.21xlO'2  3.28xlO'3  4.51xl0"6 


-9.77X10'1  1.07x1  O’4  -1.62X10'1  -3.22xl0*4  -I.66xl0'3  3.86xl0'4 
-3 .93x1  O'2  3.6x10" 5  -3.63xlO"2  -1.45xl0‘4  -3.24xl0"4  -1.29xl0"4 


JL  = 


0.44802 

-3.44841 

0.00066 


4.  Forests  and  Water 


-1.99x10" 3  -6.79x10*  -3.32xl0‘4  2.22xlO"4  3.74xl0"4  5.3xlO"4  -5.7xlO"6 
-2.75xl0‘2  7.77x10"*  2.4xl04  -7.63xl0"4  1.53xlO"2  1.52xlO"3  1.64xl0‘6 

-7.26x10'*  -4.9x1 0"7  1.14x10"*  1.9xlO"3  -4.24xl0‘3  1.61xlO"3 
6.28x10  *  1.35x10" 5  -3.5xlO"2  2.04xl0"3  -5.8xlO"4  -8.2xl0"4 


J.  = 


1.87565 

-1.05647 

0.00005 


5.  Fields  and  Water 


9.93x10* 3  -6.58x10"'  -2.98xl0'4  -2.62xlO"3  -3.99xlO"3  9.63xlO"4  -4.34xl0"6 
3.57xlO'2  -4.08x10"'  -1 .82xl0'4 -1.45xlO"3  -1.72xlO"2  -1.58xlO'3  6.17xl0"6 


-7.48X10'1  -5.76xl0‘6  8.65xlO'2  9.88xl0'5  1.05xl0'3  2.54xl0*3 
-9.05x1  O’1  -1.84x1  O'5  -1. 08x10’ 1  -4.13xl0'3  -1.5xlO'3  -2.72xl0’3 


3.90389 


a  = 


-0.26522 


-0.00047 


6.  Cities  and  Water 


A  = 


5.94xl0‘3  1.45X10’1  2.27xl0'5  -1.36xl0’3  -1.23xl0'3  -2.31xl0‘4  2.01xl0'6 
1.12xl0‘2  -9.56X10’1  -2.06x1 0*4  3.14xl0’3  -2.98xl0‘3  -8.57xlO'3  9.29xlO’7 


9.89X10’1  -2 .0x1 0‘6  -3 .5 5x1  O’2  -9.18xl0'4  5.24xl0’3  6.68xl0'4 
-2.93x1  O’1  -S.OlxlO’6  -3.51xlO*3  7.05xl0‘4  -1.29xlO'3  4.25xl0'3 


A  = 


1.31629 

0.94784 

-0.00007 


Figures  1  through  6  show  plots  of  the  training  samples  for  the  two  components  of  Y. 
Figure  1  shows  a  plot  of  the  Y  data  for  forests  and  fields.  The  X’s  represent  data  from  forest 
samples,  and  the  squares  represent  fields.  The  line  in  the  center  is  the  decision  boundary 
provided  by  the  pseudoinverse  classifier.  It  can  be  seen  that  the  data  for  forests  and  fields 
are  not  completely  separated  since  three  points  belonging  to  the  class  of  forests  are  located 
to  the  right  of  the  final  decision  boundary.  Figure  2  shows  a  plot  of  the  Y  data  for  forests 
and  cities.  These  two  classes  are  not  completely  separated  either  as  shown  by  the  location 
of  the  decision  boundary.  Figure  3  shows  a  plot  of  the  Y  data  for  cities  and  fields.  These 
two  classes  are  well  separated,  which  corresponds  to  the  large  value  of  Jj.  Figures  4,5,  and 
6  depict  the  Y  data  for  forests  and  water,  fields  and  water,  and  cities  and  water.  In  each  of 
these  three  cases  the  Y  data  are  very  well  separated  and  clearly  defined.  Once  the  decision 
boundaries  were  computed  for  each  of  the  six  possible  pairs  of  classes,  the  classifier  was 
used  to  determine  the  class  of  each  of  the  400  training  samples.  The  results  of  this  test  are 
given  in  table  2. 
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FORESTS  AND  FIELDS 


Forests  and  Fields. 


FORESTS  AND  CITIES 


.  Forests  and  Cities. 


and  Fields. 


and  Water. 


ELDS  AND  WATER 


and  Water. 


Table  2.  Results  of  Classifying  the  Original  Training  Set 


CLASS 

NUMBER 

CORRECT 

NUMBER 

WRONG 

PERCENTAGE 

CORRECT 

Forests 

94 

6 

94 

Fields 

100 

0 

100 

Cities 

99 

1 

99 

Water 

100 

0 

100 

TOTALS 

393 

7 

98.25% 

After  the  classification  techninue  was  shown  to  work  well  on  the  original  training  set, 
it  was  applied  to  the  other  300  samples  that  came  from  the  same  imagery,  but  were  not 
used  to  train  the  classifier.  The  results  of  tills  application  are  shown  in  table  3. 


Table  3.  Results  of  Classifying  the  300  Data  Samples 

NUMBER 

NUMBER 

PERCENTAGE 

CLASS 

CORRECT 

WRONG 

CORRECT 

Forests 

93 

7 

93 

Fields 

97 

3 

97 

Water 

95 

5 

95 

TOTALS 

285 

15 

95% 

It  can  be  seen  that  ihe  third-order  approach  has  worked  well  on  samples  taken  from 
the  same  imagery  that  was  used  for  training  purposes.  It  remains  to  be  seen  if  similar  results 
can  be  obtained  with  other  imagery  taken  with  the  same  radar  system. 


CONCLUSIONS 


1.  The  third-order  co-occurrence  approach  for  texture  analysis  was  capable  of  separating 
four  classes  of  terrain  features  on  radar  imagery. 

2.  The  thirteen  features  calculated  from  the  third-order  co-occurrence  array  provided  a 
98.25%  correct  classification  rate  for  the  original  training  set  and  a  95%  correct 
classification  rate  for  the  300  data  samples  that  were  taken  in  the  vicinity  of  the  training 
set. 

3.  The  excellent  correct  classification  results  obtained  in  this  study  were  also  due  to  the 
linear  feature  selection  technique  used  to  choose  features  which  were  optimized  for  showing 
class  separability. 
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Appendix  A.  Computer  Program  for  Calculating  J  j  and  the  A  Matrix 

along  with  the  Subroutine  for  Computing  the  Third-Order 
Co-occurrence. 


iSGLNC  T=00004  IS  ON  CR00023  USING  00037  BLKS  R*0000 
001  FTN4X.L 

0  0  02  C***»* ************* »***^*  ***********************************  *******C 

0003  C  THIS  PROGRAM  PERFORMS  FEATURE  SELECTION  FOR 
0004  C  MULTIDISTRIBUTIONS  BY  MAXIMIZING  THE  VALUE  OF 

0005  C  J1 .  THIS  PROGRAM  USES  THE  SPATIAL  GRAY  LEVEL  DEPENDENCE  TENSOR 
0006  C******************************************************************C 

0007  PROGRAM  SGLNC< 3 . 80 0 ) . REV . 1 0/1 3/83 

0008  DIMENSION  LUOT< 5 > . IMAGE ( 1 024 ) . X < 1 0 0 . 1 3) . IDCB< 1 44 > . IFIL2 ( 3 ) 

0009  1 . IDYL (35 .21 , A< 13 . 2) ,C < 13 . 1 3 .2) ,S1 < 1 3 . 13 > . S2( 1 3 , 13 ) .XMO ( 13) ,P (2) 

0010  2. El < 169) .LX (13) .MX< 13) .INAM<3> ,BUF<28>  . IDXL<35.2> .IFIL1 (3)  . 

0011  3  IRES <900 ) .MASK (3.3) 

0012  DATA  IANS/2HYE/ 

0013  DATA  MASK<1 . 1 1/0/ .MASK < 1 ,2>/-l/ .MASK < 1 .3)/0/ 

0014  DATA  MASK  (2.1)/-1/ . MASK ( 2 . 2  > /5/ . MASK (2.3)/  — 1/ 

0015  DATA  MASK  <3.1 )/0/.MASK(3.2)/-l/. MASK ( 3 . 3)/0/ 

0016  CALL  RMPAR (LUOT > 

0017  LL=LUOT ( 1 ) 

0018  CALL  ERLU<LU) 

0019  1  NC=2 

0020  NDATA=1 00 

0021  UR ITE ( LU . 8) 

0022  8  FORMAT < “UH I CH  LINE  PRINTER  DO  YOU  WANT  TO  USE?") 

0023  READ<LU.*)LUP 

0024  URITE(LU.IO) 

0025  10  FORMAT< "ENTER  THE  NAMES  T)F  THE  DATA  SET  FILES") 

026  READ(LU .17)  IFIL1.IFIL2 

0  027  17  FORMAT < 2 < 3A2 ) ) 

0028  URITE<LU .15) 

0029  READtLU , * >  IDLU 

0030  URITE<LU .2104) 

0031  2104  FORMAT < "HOW  MANY  SETS  DO  YOU  WANT  TO  RUN  THROUGH?") 

0032  READtLU  .*)  ISET 

0033  DO  99  1=1. ISET 

0034  DO  99  J=1 , 2 

0035  URITE<LU.98)J.J.I 

0036  98  FORMAT< "ENTER  IDX'.Il."  AND  IDY".I1."  FOR  #  ,"I3> 

0037  99  READtLU.*)  IDXLU  .  J >  . IDYL (I .  J  > 

0038  DO  61  IL=1 .ISET 

0039  URITElLUP .1330) 

0040  1330  FORMAT ( 1 X  .  "  '  A '  MATRIX  CALCULATIONS  FOR  THIRD  ORDER  CO-OCCURRENCE") 

0041  IDXl«IDXL<IL.l)  . 

0042  IDX2=IDXL<IL .2) 

0043  IDY1 =IDYL< IL . 1  ) 

0044  IDY2=IDYL<IL.2) 

0045  UR ITE <  LUP .1400)  IDX1 . IDY1 . IDX2 . IDY2 

0046  1400  FORMAT <2X . "IDX1  ■  ".13."  IDY1  »  " . I3/2X . " IDX2  *  ".13."  IDY2  =  ” 

0047  1.13) 

0048  DO  80  NJ-l.NC 

0049  15  FORMAT < "DISK  LU  NUMBER?") 

0050  IF(NJ.EO.l)  THEN 

051  WRITE(LUP ,200  >  IFIL1 

x052  CA1,L  OPENdDCB.IERR.IFILl  .0.0. -IDLU) 

0053  END  IF 

0054  IF< NJ , EQ . 2 >  THEN 

0055  WRITE<LUP ,200)  IFIL2 


non 


0056  CALL  OPENUDCB  .  IERR  ,  IFIL2, 0  .0  , -IDLU) 

0057  END  IF 

0058  200  FORMAT (IX. 3A2 ) 

0059  IFUERR.LT. 0)  THEN 

0060  WRITE(LU.2010)IERR 

061  2010  FORMAT< "OPEN  FILE  ERROR". 15) 

0062  GO  TO  999  ' 

0063  END  IF 

0064  14  ICONT=l 

0065  CALL  LABINCIDCB  .LUP) 

0066  13  J=1 

0067  DO  16  1  =  1  .8 

0068  CALL  READF ( IDCB .IERR . IMAGE (  J )  ) 

0069  16  J=J+128 

0070  IFUERR.LT. 0)GO  TO  3000 

0071  GO  TO  18 

0072  3000  ICONT=NDATA 

0073  WRITE(LU.2020)IERR 

0074  2020  FORMAT< "READ  FILE  ERROR". 15) 

0075  GO  TO  999 

0076  18  NGR AY=1 6 

0077  CALL  LOPERUMAGE. MASK  .IRES. 32) 

0078  CALL  ISCALURES. IRES. 900 .0.15) 

0  079  CALL  SGLDT  (IRES  .  SMO  .  SNE  ,  ENL  , DKN, DIN .  QJN  .  CON.  ENT  THO,. 

0  080  1 ABV  ,AQQ,A01>A10  ,  NGR  AY  ,  IDX1  ,  IDY1  ,IDX2,IDY2> 

0081  C 

0082  C======STORE  RESULTS  IN  AN  ARRAY===========================C 

0083  C 

0084  X< ICQNT . 1 )=SMO 

0085  XUCONT  ,2)=SNE 

086  XUCONT.  3>=£NL 

d 087  XUCONT  ,4)=QKN 

0088  XUCONT. 5)=0IN 

0089  XUCONT  .6)*QJN 

0090  XUCONT  ,7)=C0N 

0091  XUCONT.  8)=ENT 

0092  XUCONT. 9)=THO 

0093  XUCONT,  10>=ABV 

0094  XUCONT.  11  >=A00 

0095  X  < ICONT , 1 2 )=A01 

0096  X  < ICONT .13) =A1 0 

0097  IFUCONT-NDATA)20  ,22,22 

0098  20  ICONT  =  ICQNT  +  1 

0099  GO  TO  13 

0100  22  CALL  CLOSE(IDCB) 

0101  CALL  C0VER<NDATA,C,A,13.NJ.X.D3> 

0102  WRITEtLUP .70)D3 

0103  70  FORMATUX."  THE  INTRASET  DISTANCE* ”  .E15 . 8 ) 

0104  URITEUUP  ,140) 

0105  140  FORMATUX.  "COVARIANCE  MATRIX") 

0106  DO  142  K=1 ,13 

0107  URITE(LUP.144)<C<K.M.NJ),M=1 .13) 

Cl  08  144  FORMATUX,  13<E9. 3.  IX)) 

0109  142  CONTINUE 

0110  DO  1  000  1=1  .13 

111  DO  1  000  J=1 .13 

v*  1  1 2  1000  S2<  I  ,J)=C(  I  ,J  ,NJ> 

0113  CALL  INV<S2.13.D.LX.MX> 

0114  UR ITE ( LUP  ,  1 0 1 0 ) D 

0115  1010  FORMAT( IX . "THE  DETERMINANT  OF  THE  COVARIANCE  MATRIX*" .E15 . 8 


! 

0116 

j 

DO  1020  1-1 ,13 

■ 

0117 

DO  1020  J-1 .13 

0118 

1020  S2(I . J)*0. 0 

> 

0119 

P(l>*0,5  *  P(2>=0 .5 

i 

0120 

WRITE (LUP .  156)NJ,P(NJ) 

I 

121 

156  FORMAT <1X,"NJ=", 11 >5X, *P (NJ>«* ,F1 0 . 8) 

j 

0122 

80  CONTINUE  " 

I 

0123 

DO  90  1*1 .NC-1 

0124 

DO  90  J*I+1,NC 

1 

0125 

90  S2(I.J>*0.0 

l 

0126 

DO  72  1*1, NC-1 

M 

a 

0127 

DO  72  J*I+1  ,NC 

0128 

DO  72  M*1 . 13 

0129 

S2(  I . J )*S2( I .J)+(A(M, I )-A(M, J) )<**2 

i 

0130 

72  CONTINUE 

0131 

DO  91  1-1 .NC-1 

1 

0132 

DO  91  J-I+l . NC 

1 

0133 

S2M  . J)«SQRT(S2( I , J) ) 

0134 

91  CONTINUE 

J 

0135 

WRITE (LUP .73) 

> 

0136 

73  FORMAT (IX. “THE  INTERSET  DISTANCES*) 

0137 

DO  92  1=1 .NC-1  i 

j 

0130 

DO  92  J-I+1 .NC 

0139 

WRITE (LUP ,74) I , J ,S2( I , J)  j 

1 

0140 

74  FORMAT < IX . *I*“ ,11 ,5X . * J** ,11 ,5X , “D ( I , J )-* ,E15 . 8 ) 

J 

0141 

92  CONTINUE  ! 

3 

5 

0142 

DO  30  1-1.13 

0143 

DO  30  J-1 ,13 

0144 

S2< I . J )=0 . 0 

0145 

30  S1(I.J)*0 .0 

l 

146 

DO  32  1*1.13 

ul47 

DO  32  J-1 ,13 

1 

0140 

DO  32  K-l .NC 

1 

I 

0149 

32  S2  <  X  .  J  > -S2 ( I . J  > ♦C ( I • J , K  >  *P ( K  > 

0130 

DO  40  J-1 .13 

0151 

DO  40  M-l .13 

0152 

40  SI (J.M)-C(J.M.l >+C( J,M,2)  +  (A< J, 1 >-A( J ,2 > )*(A(M,1)-A(M,2>) 

* 

0133 

DO  41  1=1,13 

0154 

41  XMO( I )-0 . 0 

0155 

CALL  NROOT(13,Sl ,S2,XM0,EI> 

0156 

DO  163  J-1, 13 

0137 

DO  165  1-1 .13 

0158 

K=U13«(J-1  > 

0139 

S2(I .J)-EI(K) 

0160 

165  CONTINUE 

0161 

WRITE(LUP .42) 

0162 

42  FQRMAT(2X, "EIGENVALUES") 

0163 

DO  47  1-1 .13 

0164 

WRITE (LUP ,46)XM0( I) 

1 

0165 

46  FORMAT (1X.E15-8) 

! 

0166 

47  CONTINUE 

0167 

WRITE (LUP .48) 

0168 

48  F0RMAT(2X. "EIGENVECTORS") 

1 

0169 

DO  50  1*1 ,13 

0170 

WRITE (LUP .76) (S2( I , J) ,J-1 ,13) 

171 

76  FORMAT (1X,13(F8.3,1X> ) 

-172 

50  CONTINUE 

4 

0173 

64  M-2 

0174 

XJ1-0.0 

* 

* 

0175 

DO  54  1-1  .M 

21 

|  .  v ,. ’..*  ”  V/rVsV  v,V,  : 

0176 

54 

XJ1-XJ1+XMO!!) 

0177 

WRZTE1LUP  .56)XJ1 .M 

0178 

56 

FORMAT <2X . “THE  VALUE  OF  J1 - " . F 1 0 . 4 , 5X , "M« " . I 1 > 

0179 

CALL  TRMAT(S2,S1 ,13,13,0) 

0180 

MRITEtLUP .58) 

81 

58 

FOR MAT (2X. “THE  TRANSFORMATION  MATRIX  A") 

0182 

SO  60  1-1 . M  " 

0183 

WRITE (LUP .63) (SI (I .J) .J-l .13) 

0184 

63 

FORMAT OX,13(E9.3.iX) j 

0185 

60 

CONTINUE 

0186 

UR1TE(LU.3010) 

0187 

3010 

FORMAT < “DO  YOU  WANT  TO  SAVE  THE  A  MATRIX  ON  A  DISK  FILE?") 

0188 

READ ( LU .30  0 1 ) IAN 

0189 

3001 

FOR MAT (A2) 

0190 

IF! IAN . NE . IANS ) CO  TO  6100 

0191 

URXTE(LU .3002) 

0192 

3002 

FORMAT! "INPUT  NAME  FOR  DATA  FILE") 

0193 

READ  <  LU . 30  03 ) IN Ah 

0194 

3003 

FORMAT <3A2> 

0195 

3005 

FORMAT (12) 

0196 

URITE(LU .3034) 

0197 

3034 

FORMAT! "INPUT  DISK  LU  6") 

0198 

READ (LU .30  05) ICR 

0199 

CALL  CREAT! IDCB . I ERR . INAM .1,3, 0 ,-ICR  > 

0200 

IF(IERR .GE . 0>CO  TO  3100 

0201 

URITE(LU.3006)lERR 

0202 

3006 

FORMAT! "CREATE  FILE  ERROR  *".I5) 

0203 

3100 

BUFd)-FLOAT(M) 

0204 

BUF !2)*13. 

0205 

IMQ-2 

’06 

DO  3200  1-1 .M 

o207 

DO  3200  J-l, 13 

0208 

1MQ-XMQ+1 

0209 

3200 

BUF! IMS)— SI  !Z, J) 

0210 

CALL  WRITF (IDCB. IERR. BUF. 2»IMQ) 

0211 

IF ( IERR .GE . 0)GO  TO  3300 

0212 

URITE(LU . 31 07)IERR 

0213 

3107 

FORMAT! "WRITE  FILE  ERROR-". 15) 

0214 

3300 

CALL  CLOSE (IDCB. IERR) 

0215 

IF ( IERR . GE . 0 ) GO  TO  4100 

0216 

WRITE!LU. 3310) IERR 

0217 

3310 

FORMAT! "CLOSE  FILE  ERROR-". 15) 

0218 

GO  TO  999 

0219 

4100 

WRITE! LUP . 3320  > INAM . ICR 

0220 

3320 

FORMAT!/. IX. "MATRIX  A  STORED  ON  FILE  NAMED  ".3A2."  LU-",I2> 

0221 

6100 

WRITE (LUP .96) 

0222 

96 

FORMAT! "1") 

0223 

61 

CONTINUE 

0224 

999 

STOP 

0225 

END 

0226 

SUBROUTINE  COVER(KK .COV.AV.N.NJ  X.D3) 

0227 

DIMENSION  X!100 ,13) ,C(13.13>  A! 13) .COV(13 .13,2) ,AV(13 .2) .VAR (13) 

0228 

XK-KK 

0229 

DO  50  K-l .N 

0230 

A(K  >-0 .0 

"*31 

DO  45  J-l  .KK 

.232 

45 

A(K)-A!K)+X!J,K) 

0233 

50 

A(K)-A(K)/XK 

0234 

DO  55  K-l ,N 

0235 

VAR (K >-0.0 

22 


023b 

DO  54  1*1 .KK 

0237 

54 

VAR<K>-VAR<KH-<X(I.K>-A<K> >##2 

0238 

55 

CONTINUE 

0239 

DO  60  K-l ,N 

0240 

60 

VAR  <K) -VAR  <  K ) / ( XK- 1 . ) 

?41 

D3-0.0 

u242 

DO  65  K  =  1  .N 

0243 

65 

D3-D3+VAR1K) 

0244 

D3=2.*D3 

0245 

DO  1 10  K-l .N 

0246 

DO  115  1-1 .KK 

0247 

115 

X(I  ,K)«X<I .K)-A<K> 

0248 

110 

CONTINUE 

0249 

DO  120  K-l .N 

0250 

DO  125  M-l .N 

0251 

C<K.M>=0.0 

0252 

DO  130  1*1  .KK 

0253 

130 

C(K .M)*C<K.M>-*-X<I..K>*X<I.H> 

0254 

C<K ,M)«C<K.M>/XK 

0255 

125 

CONTINUE 

0256 

120 

CONTINUE 

0257 

DO  135  K«1 .N 

0258 

DO  135  M-l .N 

0259 

AV<K.NJ)*A(K) 

0260 

135 

COV<K ,H>NJ)-C<K ,M) 

0261 

RETURN 

0262 

END 

0263 

SUBROUTINE  NROOT<M.A.B,XL.X> 

0264 

DIMENSION  A ( 1 69 > ,B<169> ,XL<13> ,X(169> 

0265 

K=1 

?66 

DO  100  1*2. M 

267 

L-M*( 1-1 ) 

0268 

DO  100  I-l.J 

0269 

L=L+1 

0270 

K*x+r 

0271 

100 

B(K)-B(L) 

0272 

MV*0 

0273 

CALL  EIGEN (B.X.M. MV) 

0274 

L-0 

0275 

DO  110  J-l.M 

0276 

L-L+J 

0277 

110 

XL< J)*l . 0/SORT (ABS(B(L) ) > 

0278 

K*0 

0279 

DO  115  J-l.M 

0280 

DO  115  1-1 ,M 

0281 

K-K  +  l 

0282 

115 

B<K )*X  <K  >*XL  <  J> 

0283 

DO  120  1-1 .M 

0284 

N2-0 

0285 

DO  120  J-l.M 

0286 

Nl-M*( 1-1 ) 

0287 

L-M»< J-l )♦! 

0288 

X<L)-0 .0 

0289 

DO  120  K*1 .M 

0290 

Nl-Nl+1 

*291 

N2-N2+1 

292 

120 

X(L)-X(L)+B(N1 >»A<N2> 

0293 

L-0 

0294 

DO  130  J-l .M 

0295 

DO  130  1*1 . J 

23 


0296 

Nl-I-M 

0297 

N2*M»<J-1 ) 

0298 

L-L+l 

0299 

A(L)*0  <  0 

0300 

DO  130  K*1.M 

•01 

Nl-Nl+M 

bi02 

N2-N2+1 

0303 

130 

A<L)«A(L)+X<N1 )*B(N2> 

0304 

CALL  EIGEN <  A . X . M , MV  > 

0305 

L*0 

0306 

DO  140  1*1 .M 

0307 

L-L+I 

0308 

140 

XL(I)-A(L) 

0309 

DO  150  1-1 ,M 

0310 

N2*0 

0311 

DO  150  J-l . H 

0312 

N1*I-H 

0313 

L*M«<J-1 >+I 

0314 

A(L)*0 . 0 

0315 

DO  150  K-l.M 

0316 

N1*N1*M 

0317 

N2-N2+1 

0318 

150 

A(L>*A(L)+B(N1 )*X(N2> 

0319 

L*0 

0320 

K-0 

0321 

DO  180  J*t .M 

0322 

SUMV-0 . 0 

0323 

DO  170  1*1 .M 

0324 

L*L+1 

0325 

170 

SUMV*SUMV+ A  <L  > «A  <  L ) 

*26 

175 

SUM V*SOR T  <  SUN V ) 

.j27 

DO  180  1*1 .M 

0328 

K*K*1 

0329 

180 

X<K)*A<K  J/SUMV 

0330 

RETURN 

0331 

END 

0332 

0333 

END* 

4DSUBS  T-00004  IS  ON  CR00023  USING  00123  BLKS  R*Q000 


0682 

0683 

0684 

0683 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

0710 

0711 

0712 

0713 

0714 

0715 

0716 

0717 

0718 

0719 

0720 

0721 

0722 

0723 

0724 

0725 

0726 

0727 

0728 

0729 

0730 

0731 

0732 

0733 

0734 

0735 

0736 

0737 

0738 

0730 


SUBROUTINE  SGLDT < I AR . SMO . SNE . ENL .QKN.QIN. 0 JN . CON . ENT . THO . 
1ABV.AO0 ,401 .A10.NGRAY.IDX1.IDY1.IDX2.1DY2> 

DIMENSION  IAR (30 .30 > .COUT< 16 . 16 . 16> 

SMO*0 
SNE-0 
ENL-0 
QKN-0 
QIN«0 
OJN=0 
CON«0 
ENT=Q 
THO=0 
ABU*0 
A00*0 
A01*0 
A1 0*0 

DO  10  I-l.NGRAY 
DO  10  J-l .NCR AY 
DO  10  K-l.NGRAY 
COUT ( I .J ,K )*0 
10  CONTINUE 

IF < IDX1 . GE . 0 . AND . IDX2 . GE . 0 )GO  TO  12 
IF ( 1DX1 . LT . 0 . AND. IDX2 . GE . 0 > GO  TO  14 
IF ( IDX1 j  GE . 0 . AND . IDX2 . LT . 0  >CO  TO  16 
IFdDXl  .LT.O. AND. IDX2.LT.O>GO  TO  18 
12  NXB*1 

IF ( IDX1 . EO . 0 . AND . IDX2 . ED . 0  > NXE*30 
IFdDXl  .GE .  IDX2)NXE*30-IDX1 
IF ( IDX2 . GT . IDX 1 >NX£*30-IDX2 
GO  TO  20 
14  NXB*1-IDX1 
NXE-30-IDX2 
GO  TO  20 
16  NXB*1-IDX2 
NXE=30-IDX1 
GO  TO  20 

18  IFdDXl  .LE.  IDX2JNXB-1-IDX1 
IF (IDX2 . LT . IDX1 >NXB*1-IDX2 
NXE*30 
20  CONTINUE 

IF ( IDY1 . GE . 0 . AND , IDY2 . GE . 0  >GO  TO  22 
IFdDYl  .LT.O  .AND.  IDY2.GE.O>GO  TO  24 
IFdDYl.GE.O.  AND  .  IDY2 .  LT  .  0  )GO  TO  26 
IFdDYl  .LT.O.  AND,  IDY2.LT.  0>  GO  TO  28 
22  NYB*1 

IFdDYl  .  EG.  0  .  AND  .  IDY2  .  EO  .  0>NYE-30 
IFdDYl  .GE  .  IDY2)NYE=30-IDY1 
IF< IDY2 , GT . IDY1 >NYE*30-IDY2 
GO  TO  30 
24  NYB=1-IDY1 
NYE*30-IDY2 
GO  TO  30 
26  NYB-1-IDY2 
NYE-30-IDY1 
GO  TO  30 

28  IFdDYl  .  LE  .  IDY2)NYB*1-IDY1 
IF (IDY2.LT. IDY 1 )NYB*1 -IDY2 
NYF-30 
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0741 

DO  32  I-NY8.NYE 

0742 

DO  32  J-NXB.NXE 

0743 

Il-IAR(I ,J> 

0744 

I2-IAR< I+IDY1 .  J-MDX1 > 

0745 

I3-IAR  < I*IDY2 .  J+XDX2) 

0746 

COUT (11*1  ,12*1  , 13+1  )-C0UT< Il+l  ,12*1  .I3*1>*1  . 

0747 

32 

CONTINUE 

0748 

OC-O 

0749 

DO  40  1-1 .NCR AY 

0750 

DO  40  J-l.NCRAY 

0751 

DO  40  K-l .NCR AY 

0752 

BC-OC+COUT <  X  >J  .K) 

0753 

40 

CONTINUE 

0754 

DO  42  I-l.NGRAY 

0755 

DO  42  J-1 .NCR AY 

0756 

DO  42  K-l . NCR  AY 

0757 

XI-I 

0758 

XJ-J 

0759 

XK-K 

0760 

SMO-SMO  ♦  ( COUT  < I . J . K ) ) »»2 

0761 

SNE-SNE*COUT  < I . j . K ) / <  X I »*2*X J«#2*XK#«2  > 

0762 

ENL«ENL*<XI»«2*XJ**2*XK**2>»C0UT(I . J.K) 

0763 

CON-CON* ( <  XI-XJ  >  *»2* ( XJ-XK ) »*2*  <  XI-XK ) ««2 ) *COUT ( I . 

0764 

IE (COUT < I , J . K ) . EO . 0 . 0 ) CO  TO  44 

0765 

ENT -ENT *COUT ( I . J .K ) «AL0G1 0 ( COUT ( I . J . K ) ) 

0766 

44 

THO-THO* < COUT ( i , j , K ) > «»3 

0767 

ABU-ABV+C IABS( I-J)*IABS< J-K  >*IABS( I-K ) )*COUT < I . J ,K 

0768 

IFCIABS(I-J) • EC . 0 . AND . IABS( I-K ) . EO . 0  >400-A00+COUT ( 

0769 

IF ( I ABS ( I- J ) . EO . 0 . AND . I AB8 < I-K  > . EO . 1 > A01 -AO  1 +COUT ( 

0770 

IF ( IABS(I-J) . EO . 1 . AND . IAB8( X-K  > . EO ■ 0  > At  0-A1 0+COUT ( 

0771 

42 

CONTINUE 

0772 

DO  SO  1-1 .NCR AY 

0773 

DO  50  J-l.NCRAY 

0774 

QI-0 

0775 

QJ-0 

0776 

OK-O 

0777 

DO  52  K-t.NGRAY 

0778 

OI-flI+COUT (K . I , J) 

0779 

QJ-QJ+COUT ( IK  .  J  ) 

0780 

QK*QK*COUT(I . J.K) 

0781 

52 

CONTINUE 

0782 

QIN-0IN*QI»*2 

0783 

8JN-QJN*0v*»2 

0784 

0KN-QKN*0K»*2 

0785 

50 

CONTINUE 

0786 

SMO-SHO/QC 

0787 

SNE-SNE/OC 

0788 

ENL-ENL/QC 

0789 

ENT— ENT /QC 

0790 

THO-THO/OC 

0791 

ABV-ABU/OC 

0792 

OIN-QIN/OC 

0793 

OKN-OKN/OC 

0794 

OJN-SJN/OC 

0795 

RETURN 

0796 

END 

Appendix  B.  Computer  Program  for  Calculating  the  a.  Vector  using  the 
Pseudoinverse  Technique. 


4YAXEC  T-00004  IS  ON  CR00023  USING  00024  BLKS  R-0000 


FTN4.L 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  THIS  PROGRAM  COMPUTES  THE  TRANSFORMATION  C 

C  Y-AX  WHERE  A  IS  COMPUTED  FROM  ANOTHER  C 

C  PROGRAM  AND  IS  INPUT  HERE  C 

c  c 

C  THIS  IS  THE  PUESDO- INVERSE  METHOD  C 

C  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccnccccccccccccccc 

PROGRAM  YAX£C(3 ,300 ) 

DIMENSION  LUOT(S)  .  IMAGE  <  1  024)  .XM00.13)  . Y(2 . 1 R0 .2)  .  A<2 . 1 3) 
DIMENSION  IDCBU44) .IFILEC3) .YH<200,3) ,YT<3,200)  .YYC3.3) 

DIMENSION  LX(3>  .HX(3) , YCC3.200)  .B(200)  .AA<3>  ,l‘NA*t<3)  .DATAAC50) 

1 , INBBC3) . IDCB3(144> , IDCB2(14*) , MASK (3, 3) . IRES(?00) 

DATA  MASKU .  1 >/0/ .MASK < 1 ,2)/-l/ .MASK ( 1 ,3>/0/ 

DATA  MASK (2 >1 >/-!/  .MASK  (2 >2)/  3/  .MASK  <2  .3)-/-l/ 

DATA  MASK (3 . 1 >/0/ .MASK (3>2>/“I/ . MASK (3 .3)/0/ 

CALL  RHPAR(LUOT) 

LU-LUOT(l) 

CALL  ERLU(LU) 

Nc«2 

4  FORMAT <11 ) 

WRITE(LU.6) 

6  FORMAT ( “WHAT  IS  THE  NAME  OF  THE  FILE  CONTAINING  THE  'A'  MATRIX?*) 
READ(LU.22)INAM 
WRXTE(LU, 15) 

READ <LU,») ICR 

CALL  OPEN < 1DCB2 . IERR . XNAH ,0.0. ICR) 

IF ( IERR . CE . 0 )  GOTO  17 

WRITE(LU,2010)IERR 

STOP 

17  CALL  READF ( IDCB2 . IERR , DATAA > 

IF ( IERR . GE . 0  >  GOTO  70 
WRITE<LU. 2020) IERR 
STOP 

70  IMQ-2 

M-INT<DATAA<1)) 

DO  3001  X-l.M 
DO  3001  J-1 .13 
IMQ-XMQ+1 

3001  A<r.J)-DATAA(XMQ) 

WRITE (6. 501 ) 

301  FORMAT (SY. . 'SOLUTION  FOR  THE  'A'  VECTOR -PSEUDO  INVERSE  TECHNIQUE - 

1 NORTH  CAROLINA  IMAGERY".//) 

WRITE <6, 79) 

79  FORMAT < 1 X . "THE  TRANSFORMATION  MATRIX  A") 

DO  78  1-1 .M 

WRITE<6,77) < A< I , J) .J-1 .13) 

77  FOR  MAT (1X.13(E9,3.iX)) 

78  CONTINUE 

DO  34  NJ-1 .2 
DO  34  1-1.100 
DO  34  K-1.2 
34  Y<NJ.I.K)-0 
DO  80  NJ-1 .NC 
WRITE <LU. 16) 

16  FORMAT ("ENTER  THE  NUMBER  OF  IP  AGE  SAMPLES  TO  BE  ANALYZED  LE.100") 
READtLU. 18>NDATA 


0059 

18 

FORMAT  < 13) 

0060 

URITEtLU . 20 ) 

0061 

20 

FORMATt "ENTER  THE  FILE  NAME  FOR  THE  DATA  SET") 

0062 

READtLU .22) IF1LE 

0063 

22 

FORMAT <3A2> 

64 

URITEtLU. 15) 

w  «#  65 

15 

FORMAT ("DISK  LU  NUMBER?") 

0066 

READtLU .21 01 ) IDLU 

0067 

2101 

FORMAT  < 12) 

0068 

IF(NJ.CT.l)  GOTO  2204 

0069 

UR ITEtLU . 21 00 ) 

0070 

2100 

FORMAT ( "ENTER  VALUES  FOR  IDX1 . IDY1 . IDX2 . IDY2* ) 

0071 

READ ( LU . * ) IDX1 ,IDY1 ,  IPX2.IDY2 

0072 

2204 

WRITE( 6 . 1400 ) IDX1 .IDY1 .IDX2.IDY2 

0073 

1400 

FORMAT  <2X . "IDX1-" .  12 .5X . "IDY1-" . 12 .5X . "IDX2-" . 12 .5X . "IDY2-" .12) 

0074 

WRITE (6 .200 ) IFILE 

0075 

200 

FORMAT <ix.3A2) 

0076 

CALL  OPENt IDCB  .  IERR . IFILE . 0 . 0 . -IDLU ) 

0077 

IF ( IERR . LT . 0 ) CO  TO  2000 

0078 

GO  TO  24 

0  079 

2000 

URITE(LU .201 0)1 ERR 

0080 

2010 

FORMATt "OPEN  FILE  ERROR". 15) 

0081 

GO  TO  999 

0082 

24 

ICONT-1 

0083 

CALL  LABINtlDCB ,6) 

0084 

13 

J«1 

0085 

DO  19  I»I .8 

0086 

CALL  READFt IDCB. IERR. IMAGE(J)) 

0087 

19 

J-J+12B 

0088 

IF ( IERR . LT . 0 ) GO  TO  3000 

•“*99 

GO  TO  26 

>0 

3000 

ICONT-NDATA 

0091 

URITEtLU. 2020)IERR 

0092 

2020 

FORMATt "READ  FILE  ERROR". 15) 

0093 

GO  TO  999 

0094 

26 

NCR AY-1 6 

0095 

CALL  LOPERt IMAGE. MASK  IRES. 32) 

0096 

CALL  ISCALtIRES, IRES. 900, 0.15) 

0097 

CALL  SGLDT ( IRES . SMO . SNE . ENL . OKN , GIN . Q JN . CON . ENT . THO . 

0098 

1ABU.A0O.A01 .A10.NGRAY.IDX1 .IDY1 , IDX2 , IDY2) 

0099 

X( ICONT . 1 )-SMQ 

0100 

XtlCONT ,2) -SNE 

0101 

X( ICONT. 3) -ENL 

0102 

XtlCONT . 4) -OKN 

0103 

X( ICONT. 5)-0IN 

0104 

XtlCONT . 6)-0JN 

0105 

XtlCONT. 7) -CON 

0106 

XtlCONT .8) -ENT 

0107 

X( ICONT, 9) -THO 

0108 

XtlCONT ,10)-ABV 

0109 

X(  ICONT.  1D-A00 

0110 

XtlCONT .12)-A01 

0111 

X( ICONT. 13>-A10 

0112 

IFt ICONT-NDATA) 28 .30 . 30 

0113 

28 

ICONT-ICONT+1 

0*  »4 

GO  TO  13 

d  .5 

30 

URITEtLU. 32)NJ. IFILE 

0116 

32 

FORMATt "NJ-". 11 , 3X.3A2) 

0117 

DO  36  1-1 .NDATA 

0118 

DO  36  K-l . M 

28 


0119  DO  36  MK=1 .13 

0120  Y  <NJ.  I  .K)-Y!NJ.I.K  >*A!K  .MK)*X!I.MK> 

0121  36  CONTINUE. 

0122  80  CONTINUE 

0123  DO  82  1*1  .NDATA 

24  YH! 1 . 1 >*1 

ll  125  YH(I.2)*Y<l.I.l) 

0126  YH( I , 3 ) *Y  < 1 . 1 . 2) 

0127  82  CONTINUE 

0128  NT*2»NDATA 

0129  N1*NDATA+1 

0130  DO  84  I*N1 .NT 

0131  IQ-I-NDATA 

0132  YH(I.l)  —  1 

0133  YH(I «2>*-Y(2.IO  .1) 

0134  YH< 1 .3  >*-Y!2 .10.2) 

0135  84  CONTINUE 

0136  CALL  TRMAT! YH .YT.200.3.0) 

0137  CALL  PRHAT  ! YT , YH  > YY , 3 , 200 ,3) 

0138  CALL  INV! YY . 3 .D , LX . MX) 

0139  WR1TE!6,86)D 

0140  86  FORMAT! IX. “THE  DETERMINANT  OF  YY-".E15.8> 

0141  CALL  PRMAT  ! YY  f  YT. YC . 3.3 . 20  0 ) 

0142  DO  88  1-1.200 

0143  B! I>“1 

0144  88  CONTINUE 

0145  DO  90  1*1 .3 

0146  AA( I )*0 

0147  90  CONTINUE 

0148  DO  92  1*1 .3 

49  DO  92  J*1 .200 

ii  1 50  AA ! I  )*AA(  I  )>YC  !I.J)*B!J) 

0151  92  CONTINUE 

0152  WRITE!6 . 94) 

0153  94  FORMAT < i X . "THE  WEIGHT  VECTOR  A  IS  GIVEN  BELOW" ) 

0154  DO  96  1*1.3 

0155  WR I TE ( 6 . 98 ) AA ! X ) 

0156  98  FORMAT! 1X.E1S. 8) 

0157  96  CONTINUE 

0158  WRITE!LU .301 0 ) 

0159  3010  FORMAT! “WOULD  YOU  LIKE  TO  SAVE  THE  'A'  VECTOR  ON  A  DISK  FILE? " ) 

0160  .  READ !LU.14) IANS 

0161  14  FORMAT ! A2) 

0162  IF ! IANS . NE . 2H YE  >  GOTO  999 

0163  WRITE ! LU .3002 ) 

0164  3002  FORMAT! "INPUT  A  NAME  FOR  THE  FILE") 

0165  READ!LU.22)INBB 

0166  WRITE !LU .15) 

0167  READ !LU . *) ICR 

0168  CALL  CREAT ! IDCB3 . IERR . INBB .1 .3 . 0 .-ICR > 

0169  IF ! IERR . GE . 0  >  GOTO  1165 

0170  URIT£(LU .3006) IERR 

0171  3006  FORMAT! “CREATE  FILE  ERROR  4#  ".14) 

0172  STOP 

0173  1165  CALL  WRITF ! IDCB3 . IERR , AA . 6 > 

74  IF(IERR.GE.O)  GOTO  3310 

w  .75  UR ITE (LU .3107)  IERR 

0176  3107  FORMAT! "WRITE  FILE  ERROR  *:  ".14) 

0177  3310  CALL  CLOSE !IDCB3> 

0178  WRITE! 6 .3320 ) INBB . ICR 

0170  3320  F0P"AT.J-.  "The  '  A  •  VECTOR  IS  STORED  ON  "Il_E  NAMED  "  THE  “  _  - 

01 80  HE  ")"> 

0181  999  STQ>. 

0*82  END 

<*183  END* 
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Appendix  C.  Computer  Program  for  Classifying  Samples. 


6CLSEC  T-00004  IS  ON  CR00023  USING  00020  BLKS  R-0000 


0001 
0002 
"003 
004 
0005 
0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 
0025 
0026 
0027 
0028 
029 
0030 
0031 
0032 
0033 
0034 
0035 
0036 
0037 
0038 
0  039 


FTN4.L 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  THIS  PROGRAM  CLASSIFIES  SAMPLES  BASED  ON  THE 
C  MINIMUM  SQUARED  ERROR  PSEUDINUERSE  TECHNIQUE 

ccccccccccccccccccccccceccccccccccccccccccccccccccccccccccccccccccccccc 

PROCRAM  CLSECC3 .300  ) 

DIMENSION  LUOT(5) .IMAGE (1024) ,X< 13) .XQ<3> .I8TG(2> .INAM<9.6> . 
1IDCBC144) , IDCF (144) , IFILE(3) .IFIL (3) ,ATEN<6 j2,13)  .U(6.3> 

2 . XT  <  6  > .DATAAC50 ) . ISP AC< 4 ,6) . IRES ( 900 > .MASK (3.3) 

DATA  I N AM/2HF 0 , 2HR E , 2HS  T . 2HS  , 2HAN . 2HD  . 2HFI . 2HEL , 2HDS 
1 .2HF0.2HRE.2HST.2H8  .2HAN.2HD  . 2HCI ,2HTI .2HES 
2.2HCI .2HTI .2HES.2H  A.2HND.2H  F , 2HIE . 2HLD  ,2HS 
3 .2HF0 .2HRE . 2HST.2HS  .2HAN.2HD  . 2HUA . 2HTE . 2HR 
4 .2HFI . 2HEL . 2HDS.2H  A,2HND,2H  U.2HAT .2HER . 2H 
4.2HCI.2HTI.2HES.2H  A.2HND.2H  U . 2HAT , 2HER . 2H  / 

DATA  1STC/2HA  . 2HW  / 

DATA  MASK (1.1 )/0/ .MASK < 1 .2>/-l/,MASK <1 .3)/0/ 

DATA  MASK (2, 1 )/-l / .MASK <2,2)/5/  J1ASK (2,3)/— 1/ 

DATA  MASK <3.1 )/0/ .MA8K (3 .2 )/~l/ .MASK (3 ,31/0/ 

CALL  RMPAR(LUOT) 

LU-LUOT(l) 

CALL  ERLU(LU) 

ISPACd . 1 >-l  *  ISPAC(2.1)-0  0  ISP AC (3.1 )-2  •  ISPAC(4.1>«0 

ISPACC1 .2)*2  «  ISP AC <2.25*0  t  IBPAC(3.2)*2  •  ISPAC(4.2>*1 

ISPAC(1.3>-1  «  ISP AC(2 .35*0  O  ISPAC(3,3)-2  *  ISPACC4.35-0 

ISPAC(1 ,4)»1  «  I8PAC<2 .41*0  O  ISPAC<3.4>*2  •  I8PAC<4.4)*0 

ISPACd.5)«l  «  ISPAC(2 .5)*0  t  ISPAC<3,5)*2  *  ISPAC(4.5)*0 

ISPACC1 ,6)*1  *  1SPAC(2,6)“0  •  ISPAC(3.6)*2  *  ISPAC(4.6)-0 

DO  110  J-1.2 
DO  110  1-1,6 

WRITE ( LU . 1 )  ISTG( J) . ( INAM(K , I > . K-l .9) 

1  FORMAT< "ENTER  THE  FILE  NAME  CONTAINING  THE  *.A2. "MATRIX  FOR  ",9A2> 
1010  READ(LU.16)  IFIL 
WRITE(LU.18> 

READ(LU.v) ICR 

CALL  OPEN ( IDCF . I ERR . IFIL , 0 , 0 , ICR  > 

IF(IERR.GE.O)  GOTO  40 
UR ITE ( LU .2010)  I ERR 
STOP 


0040 

40 

CALL  READF(IDCF.XERR.DATAA> 

0041 

IF ( IERR . GE . 0  >  GOTO  60 

0042 

UR ITE (LU, 2020  > IERR 

0043 

STOP 

0044 

60 

IFU.EQ.2)  GOTO  66 

0045 

M«INT(DATAA(1 )) 

0046 

IMQ-2 

0047 

DO  33  IU-1  .M 

0048 

DO  33  JU-1 ,13 

0049 

IMQ-IMQ+1 

0050 

33 

ATEN ( I . IU . JU ) -DATAA  < IMO ) 

0051 

GOTO  119 

0052 

66 

DO  35  IU-1. 3 

0053 

35 

W(1  .IU)-DATAA(IU) 

054 

119 

CALL  CLOSE (IDCF) 

0055 

110 

CONTINUE 

0  056 

5 

WRITE(LU.IO) 

0057 

10 

FORMAT < "ENTER  THE  NUMBER  OF  IMAGE  SAMPLES  TO  BE  ANALYZED. LE.  100") 

0058 

READ(LU . 12)NDATA 

30 


o  n 


0039 

0060 

0061 

0062 

0063 

64 

u065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0073 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0083 

0086 

0087 

0088 

189 

*090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0103 

0106 

0107 

0108 

0109 

0110 

0111 

0412 

0113 

114 

115 
0116 
0117 
0118 


12  FORMAT < 13) 

URITE(LU .14) 

14  FORMAT< "ENTER  THE  NAME  OF  THE  DATA  SET") 

READ (LU , 16) 1FILE 
16  FORMAT (3A2) 

URITE(LU. 18) 

18  FORMAT ("DISK  LU  NUMBER?") 

READ ( LU .2101) IDLU 
2101  FORMAT  < 12) 

20  FORMAT ( IX. 3A2) 

CALL  OPEN ( IDCB , IERR , IFILE ,0.0, -IDLU ) 

IF (IERR  .  GE.O)GO  TO  22 
2000  WRITE (LU .201 0 ) IERR 
2010  FORMAT ( "OPEN  FILE  ERROR*. I5> 

60  TO  999 
22  IC0NT-1 

WRITE(6.20 ) IFILE 
CALL  LABIN1IDCB  >6) 

13  J»1 

DO  24  1-1,8 

CALL  READF(IDCB,IERR..IMAGE(J>> 

24  J-J+128 

IF( IERR. 6E. 0)60  TO  26 
3000  ICONT-NDATA 

WRITE(LU. 2020) IERR 
2020  FORMATCREAD  FILE  ERROR ".15) 

60  TO  999 
26  N6RAY-16 

CALL  LOPER( IMAGE. MASK. IRES. 32) 

CALL  I8CAL(IRES> IRES, 900 ,0,15) 

DO  300  1-1,6 
IDXl-ISPAC(l.I) 

IDY1-I8PAC(2,I) 

IDX2-I8PAC(3,D 

CALL  SCLDTdRES  ,X(  1 > ,X(2) ,X(3> .X(4> .X(5) .X(6> ,X(7>  fX(8> ,X(9> , 
1X00), X(11)  .X(  12)  ,X(  13)  .N6RAY.IDX1  .IDY1  .IDX2.IDY2) 

34  Yl-0 
Y2-0 

DO  36  7-1.13 
Yl«Y1*ATEN(I..l  ,J)»X(J) 

Y2-Y2«-ATEN(I .2, J)«X(7> 

36  CONTINUE 
XO( 1 )-l 
XO(2)-Y1 
XO(3>«Y2 
XT (l)-0 
DO  38  J-l .3 

38  XTd  )-XT  ( I  )+W(I .  J)*XO(  J) 

300  CONTINUE 
KONTl-O 
KONT2-0 
KONT3-0 
KONT4-0 

IF(XT( 1 ) . GT . 0 i 0  >K0NTl -KONT1+1 
IF(XT(1 ) .LT.O . 0)KONT2»KONT2«T 
IF ( XT (2 > .  GT .  0 . 0  )K0NT1  -KQNTIO 
IF (XT (2) .LT , 0 . 0  )K0NT3-KQNT3+1 
IF(XT(3> . GT . 0 . 0  >K0NT3-K0NT3+1 
IF (XT (3) .LT.O. 0)KONT2«KQNT2*1 


IF(XT(4> ,GT.O.O)KONTI-KONT1+1 
IF  (XT (4 )  .LT.O  . 0)KONT4-KONT4+1 
IF(XT<5> . GT . 0 . 0 )K0NT2*K0NT2+1 
IF < XT (5) . LT . 0 . 0 )KONT  4*K  QNT4+1 
IF(XT(6) ,GT.O.O)KQNT3-KONT3-M 
IF (XT (6) .LT.O .0>KONT4«KONT4*1 

IF  (KONT1  .  CT .  K0NT2“?AND .  KQNT1  . GT . KDNT3 .  AND  .  KQNT1  .  GT .  KQNT4  )GO  TO  72 
IF(K0NT2 . CT . KONTI .AND . KQNT2 . GT . K0NT3 . AND . KQNT2 . CT . K0NT4) GO  TO  74 
IF < K0NT3 . CT . XONT1 . AND . K0NT3 . CT . K0NT2 . AND . KONT3 . CT. KONT4)CO  TO  76 
IF (K0NT4 . CT . KONTt . AND . K0NT4 . CT . KQNT2 . AND . K0NT4 . CT . K0NT3 >GO  TO  78 
WRITE(6 .302) 

302  FORMAT ( IX . "UNDETERMINED” ) 

GO  TO  90 
72  URIT£(6 .80  >  . 

80  FORMAT (IX. "FOREST “) 

GO  TO  90 
74  WRITE (6. 82) 

82  FORMAT (IX. "FIELDS") 

GO  TO  90 
76  URITE(6 .84) 

84  FORMAT (IX. "CITY") 

CO  TO  90 
78  WRITE(6 . 86) 

86  FORMAT (IX. "WATER") 

90  IF ( ICONT-NDATA) 92 .94 . 94 
92  ICONT-ICONT+1 
GO  TO  13 

94  CALL  CLOSE (IDCB) 

WRITE(LU .995) 

995  FORMAT ( "DO  YOU  WANT  TO  GO  THROUGH  AGAIN?*) 

READ (LU. 994) I AS 
994  FORMAT (A2> 

IF ( IAS . EG . 2HYE)  GOTO  5 
999  STOP 
END 
END* 


WWP 


mm 


